Image taken from http://www.martinos.org/neurorecovery/technology.htm
cort_thickness in extrantsr packagelabel roi 1 X1002 “left caudal anterior cingulate” 2 X1003 “left caudal middle frontal” 3 X1005 “left cuneus” 4 X1008 “left inferior parietal” 5 X1017 “left paracentral” 6 X1022 “left postcentral” 7 X1023 “left posterior cingulate” 8 X1024 “left precentral” 9 X1025 “left precuneus” 10 X1027 “left rostral middle frontal” 11 X1028 “left superior frontal” 12 X1029 “left superior parietal” 13 X1030 “left superior temporal” 14 X1031 “left supramarginal” 15 X1035 “left insula”
source('utils.R')
source('combat.R')
modelData = read.csv('modelData.csv')
head(modelData)
subject age sex dx scanner 1 002_S_0413 76.4329 F Normal 3_PHILIPS_Intera_2 2 002_S_0559 79.4658 M Normal 3_PHILIPS_Intera_2 3 002_S_0729 65.2986 F MCI 3_PHILIPS_Intera_2 4 002_S_0816 70.9534 M AD 3_PHILIPS_Intera_2 5 002_S_0954 69.5041 F MCI 3_PHILIPS_Intera_2 6 002_S_1018 70.8055 F AD 3_PHILIPS_Intera_2
mod = model.matrix(~age+factor(sex)+factor(dx), data=modelData)
ctData = read.csv('imageData.csv')
head(ctData)[,1:10]
subject X1002 X1003 X1005 X1008 X1017 X1022
1 002_S_0413 2.349515 1.957340 1.6771022 2.929418 1.893073 1.733991
2 002_S_0559 2.814481 1.768518 0.9173002 2.297948 1.612640 2.071636
3 002_S_0729 2.788202 2.108902 1.6343228 3.154604 2.219242 1.984391
4 002_S_0816 2.753477 2.168249 1.7955870 2.880065 1.973897 2.042263
5 002_S_0954 2.523273 1.664635 1.0189855 2.077749 1.236037 1.468252
6 002_S_1018 2.596688 2.216399 1.9206076 2.424231 1.744081 1.800574
X1023 X1024 X1025
1 2.365209 1.968747 2.474627
2 2.606214 1.838712 2.345573
3 2.578613 1.782150 2.055732
4 2.188748 1.990839 2.992091
5 1.329828 1.468782 1.623318
6 2.165901 1.947617 2.611285
img = t(ctData[,-1]) head(img)[,1:10]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
X1002 2.349515 2.8144805 2.788202 2.753477 2.523273 2.596688 1.833424
X1003 1.957340 1.7685180 2.108902 2.168249 1.664635 2.216399 1.696198
X1005 1.677102 0.9173002 1.634323 1.795587 1.018986 1.920608 1.228951
X1008 2.929418 2.2979484 3.154604 2.880065 2.077749 2.424231 2.712218
X1017 1.893073 1.6126404 2.219242 1.973897 1.236037 1.744081 1.760396
X1022 1.733991 2.0716363 1.984391 2.042263 1.468252 1.800574 1.712189
[,8] [,9] [,10]
X1002 2.960386 3.191757 1.6549346
X1003 2.220978 2.294089 1.1845660
X1005 1.638375 1.760091 0.6901224
X1008 2.727885 2.280574 1.2541459
X1017 2.105533 1.890207 1.0815922
X1022 2.095173 1.480126 0.8755437
harmonized = combat(dat=img, batch=modelData$scanner, mod=mod)
[combat] Performing ComBat with empirical Bayes [combat] Found 14 batches [combat] Adjusting for 4 covariate(s) or covariate level(s) [combat] Standardizing Data across features [combat] Fitting L/S model and finding priors [combat] Finding parametric adjustments [combat] Adjusting the Data
combat() function returns a list.head(harmonized$dat.combat)[,1:10]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
X1002 2.075969 2.5403388 2.514211 2.479313 2.2494910 2.322785 1.5601570
X1003 1.604147 1.3928650 1.774367 1.829854 1.2947712 1.890831 1.3247258
X1005 1.251037 0.6043039 1.215196 1.363586 0.6807333 1.459647 0.8726338
X1008 2.242658 1.6391016 2.457074 2.191413 1.4269869 1.755756 2.0304396
X1017 1.414838 1.1241720 1.749303 1.497199 0.7381837 1.263199 1.2756135
X1022 1.172603 1.5470296 1.459725 1.521212 0.8741849 1.260857 1.1429860
[,8] [,9] [,10]
X1002 2.686305 3.134451 1.7216344
X1003 1.886106 2.388280 1.3679891
X1005 1.219925 1.922309 0.8532241
X1008 2.053663 2.445597 1.4071359
X1017 1.632919 2.067534 1.2650911
X1022 1.579223 1.739706 1.1057018
harmonized$gamma.hat versus harmonized$gamma.starharmonized$delta.hat versus harmonized$delta.star